LiDAR (Light detection and ranging) is a data acquisition technology that uses light emitted from a laser pulse (up to 400000 pulses of light per second) to measure the travel time of the pulse from the source to the target and back (returns). Using this approach the instrument collects three dimensional data in large volumes, high density and at unprecedented accuracy (Figure 1). The resolution of the data collected depends on the pulse frequency, that is, resolution increases with increasing pulse frequency.
LiDAR data is provided in LAS format and classified into ground and non-ground (vegetation) returns. This type of classification is usually performed by the data provider.
Canopy density map -. Forest canopy density is the ratio of the points classified as vegetation to the point classified as ground. This variable, together with canopy height, are commonly used for biomass estimation and the assessment of vegetation coverage.
Digital Terrain Model (DTM) -. The DTM represents the bare ground, it can be calculated from the average value of points that have been classified as ground in each cell.
The package rgrass7 provides the interface between GRASS GIS 7 and the R software that allows the use of all the GRASS GIS commands from the R command line. To load the rgrass7 package into your R session, it is necessary to first install the package using the command install.package(“rgrass7”) and then:
library(rgrass7)
## Loading required package: XML
## GRASS GIS interface loaded with GRASS version: GRASS 7.0.4svn (2016)
## and location: DIARS
It is recommended to first create a folder where the datasets (LAS files, download a sample dataset from HERE) are located and the outputs and results will be saved. To set the working directory in R to the folder created:
setwd("/home/garzonc/Desktop/DIARS/")
Now you need to create a directory that will contain the GRASS GIS database for the working session:
# dir.create('grassdata')
To initiate GRASS into your R session you need to set the gisBase argument to GRASS binaries and libraries directory path, this should be something like C:\Program Files\GRASS GIS 7.00beta4 if your are using Windows (check HERE for other operating systems). This will create a new GRASS location called Sylt_island inside the GRASS GISDBASE.
myGRASS <- "/usr/local/grass-7.0.4svn"
myGISDbase <- "/home/garzonc/Documents/grassData"
myLocation <- "DIARS"
myMapset <- "PERMANENT"
initGRASS(myGRASS, home = tempdir(), gisDbase = myGISDbase,
location = myLocation, mapset = myMapset,
override = TRUE)
## gisdbase /home/garzonc/Documents/grassData
## location DIARS
## mapset PERMANENT
## rows 1440
## columns 1440
## north 6075220
## south 6074500
## west 453999.9
## east 454719.9
## nsres 0.5
## ewres 0.5
## projection +proj=utm +no_defs +zone=32
## +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000
## +to_meter=1
And now set the projection system, for this example we need WGS84 UTM zone 32 (epgs: 32632). You can search for the epgs code of your projection system here:
execGRASS("g.proj", flags = c("c"), parameters = list(epsg = 32632))
## Projection information updated
To set the geographic region in GRASS including the extent of the bounding box (i.e. the region where all the calculations will be performed) and the resolution:
execGRASS("g.region", parameters = list(n = "6075220.40",
s = "6074500.40", e = "454719.89", w = "453999.89",
res = "0.5"))
To make sure that the GRASS location is correct you can check the location metadata:
gmeta()
## gisdbase /home/garzonc/Documents/grassData
## location DIARS
## mapset PERMANENT
## rows 1440
## columns 1440
## north 6075220
## south 6074500
## west 453999.9
## east 454719.9
## nsres 0.5
## ewres 0.5
## projection +proj=utm +no_defs +zone=32
## +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000
## +to_meter=1
Once the GRASS in R environment is set you can import the LiDAR data. The following for loop will create 4 vector point files from the binary LAS files using the v.in.lidar GRASS module. Note that this will take a couple of hours to run depending on your computer configuration because this files contain altogether 16568266 points.
# for (i in
# c('GLOBAL-H_454000_6074500.las','GLOBAL-H_454000_6075000.las','GLOBAL-H_454500_6074500.las','GLOBAL-H_454500_6075000.las'))
# {execGRASS('v.in.lidar',
# flags=c('r','o','overwrite','quiet'),
# parameters=list(input=i,output=substr(i,8,23)))}
Now combine the 4 imported vector map layers into one single vector map layer (this step is also a time consuming task):
# execGRASS('v.patch', flags=c('e',
# 'overwrite','quiet'),
# parameters=list(input=c('H_454000_6074500','H_454000_6075000','H_454500_6074500','H_454500_6075000'),output='las_file'))
Select the points classified as “vegetation” and create a new vector map:
# execGRASS('v.extract',flags=c('overwrite','quiet'),
# parameters=list(input='las_file',where='class
# IN ('Low Vegetation','High
# Vegetation','Medium
# Vegetation')',output='las_file_vegetation'))
Select the points classified as “ground” and create a new vector map:
# execGRASS('v.extract',flags=c('overwrite','quiet'),
# parameters=list(input='las_file',where='class
# IN
# ('Ground')',output='las_file_ground'))
Select the “non-vegetation” points and create a new vector map:
# execGRASS('v.extract',flags=c('overwrite','quiet'),
# parameters=list(input='las_file',where='class
# IN ('Unclassified','Ground','Building'
# ,'Water')',output='las_file_non_vegetation'))
Convert each of the 3 GRASS binary vector map containing the vegetation, ground and non-vegetation points into a GRASS ASCII vector map:
# execGRASS('v.out.ascii',flags='overwrite',parameters=list(input='las_file_vegetation',output='las_file_vegetation',format='point'))
# execGRASS('v.out.ascii',flags='overwrite',parameters=list(input='las_file_non_vegetation',output='las_file_non_vegetation',format='point'))
# execGRASS('v.out.ascii',flags='overwrite',parameters=list(input='las_file_ground',output='las_file_ground',format='point'))
To visualise a subset of the point cloud using the RGL package you need to first install the rgl package using the install.package("rgl") command:
library(rgl)
ground <- read.table("las_file_ground", sep = "|",
dec = ".")
ground <- ground[which(ground[, 1] < 454100 &
ground[, 2] > 6074800 & ground[, 2] <
6074900), ]
plot3d(ground[, 1], ground[, 2], ground[,
3], col = "coral", aspect = c(1, 1, 0.1),
xlab = "", ylab = "", zlab = "", box = F,
axes = F)
vegetation <- read.table("las_file_vegetation",
sep = "|", dec = ".")
vegetation <- vegetation[which(vegetation[,
1] < 454100 & vegetation[, 2] > 6074800 &
vegetation[, 2] < 6074900), ]
plot3d(vegetation[, 1], vegetation[, 2],
vegetation[, 3], col = "green", aspect = c(1,
1, 0.1), add = T)